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Abstract 

For some anisotropic wave models, the PML (perfectly matched layer) 
method of open boundaries can become polynomially or exponentially un- 
stable in time. In this work we present a new method of open boundaries, 
the phase space filter, which is stable for all wave equations. 

Outgoing waves can be characterized as waves located near the bound- 
ary of the computational domain with group velocities pointing outward. 
The phase space filtering algorithm consists of applying a filter to the 
solution that removes outgoing waves only. 

The method presented here is a simplified version of the original phase 
space filter, originally described in [23] for the Schrodinger equation. We 
apply this method to anisotropic wave models for which the PML is un- 
stable, namely the Euler equations (linearized about a constant jet flow) 
and Maxwell's equations in an anisotropic medium. Stability of the phase 
space filter is proved. 

1 Introduction 

We consider the solution of linear wave equations, of the form: 

ut{x,t) — T-Lu{x,t) (1-1) 

Here, u : M" (or — > C") and 7i is a skew-adjoint linear differen- 

tial operator. This is a prototypical example of a linear wave equation. Such 
equations admit linear waves of the form e*('=-=^-"j('=)*)4(fc) with dj{k) the j'th 
eigenvector of Ti and (^j{k) the j-th eigenvalue of Ti. in the frequency domain. 

A great many systems of practical importance can be recast in such a form, 
including the wave equation, Schrodingers equation, the linearized Euler equa- 
tions. Maxwell's equations and others. Very often, the equations are too compli- 
cated to solve exactly, so numerical approximations must be used. If we are only 
interested in the solution of u{x, t) in some finite region B C M^, it simplifies 
the computations to solve (|l.ip only on this region. Of course, steps must be 
taken to prevent spurious reflection from the artificial boundary. 
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Exact non-reflecting boundary conditions (NRBCs) can be constructed for 
many wave equations (TJ [U [8l [9l [TOl [131 [15]. Such methods are accurate, but 
non-local in time and space, and algorithms are sensitive to the shape of the 
boundary and the particular interior solver. They also prevent the use of Fourier 
spectral methods for solving the interior problem, since these methods have 
periodic or Dirichlet boundaries built in. 

The Perfectly Matched Layer (PML) is another approach to the problem 
[31[S]. The PML consists of the complex change of coordinates, x i-^ x + ia{x) 
(for a suitably chosen (j{x)) where a{x) = Q lov x € B (this technique was used 
in scattering theory for a long time [TTJ [ID])- Inside B the dynamics of (jl.ip 
are unchanged, while on B^ plane waves are turned into exponentially decaying 
modes. The domain of computation is taken to be a region Be D B which is 
B surrounded by an absorbing layer. The layer dissipates outgoing waves and 
introduces minimal spurious reflections. Unlike NRBCs, the PML is local in 
space and time and compatible with fast spectral methods. For this reason, the 
PML is the method of choice in modern wave propagation. 

However, the PML has a fatal flaw which limits its use in some cases. For 
certain types of anisotropic wave, the PML can become exponentially unstable 
in time regardless of the particular method used (the instability exists on the 
level of the PDF) . This was first noticed in [TJ] for the Euler equations linearized 
about a jet flow (see also [11[31[IS]). In [3], a simple geometrical criterioifl was 
provided for a PML to be unstable. Consider a boundary at xi = L and a 
dispersion relation cu{k): if fci > while the group velocity dkiU}{k) < for 
some k, then the PML is exponentially unstable. 

In [22l [23] a new approach to the problem of open boundaries was intro- 
duced, the Time Dependent Phase Space Filter (TDPSF). Phase space analysis 
is a major tool in modern scattering theory [ITl [HI [19] . The key idea in this 
approach is that certain regions of quantum phase space (the set of points 
{{x,k) G M.^ X M^}, where x represents a position and k a spatial frequency) 
consist solely of outgoing waves, whereas other regions have more complicated 
interactions. The philosophy of the TDPSF is to identify parts of the solution 
u{x, t) localized in the outgoing regions of phase space, and filter them from the 
solution before they reach the boundary. After these waves are filtered, u{x, t) is 
not approaching the boundary, and therefore boundary conditions don't matter. 

In [211 [m US] , this approach is used to construct open boundaries for the 
Schrodinger equation. The phase space projections there are based on the Gaus- 
sian windowed Fourier transform [7]. In this paper, we extend the work of 
[211 [22l [23] to more general wave equations including the linearized Euler equa- 
tion. We also simplify the method significantly, replacing the windowed Fourier 
transform by standard phase space projections of the form x{x)P{k)x{x) (this 
approach was also taken in [H]). We begin by briefiy reviewing the dynamics 
of linear waves (Section II. ip . In Section [2] we present the phase space filtering 
algorithm for general systems of the form (|l.ip . In Section [3] we show some nu- 

^ There are some other conditions necessary for the proof to apply, but it is convincingly 
argued that this is the fundamental criteria. 
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merical examples. In Section HJ we briefly discuss a method of filtering outgoing 
waves with frequencies too smaU to resolve on an absorbing boundary layer of 
reasonable width. 



1.1 Dynamics of Waves 

For concreteness and to pin down our notation, we review a few facts about wave 
propagation. Recall that we defined u)j{k) as the fc'th branch of the dispersion 
relation, and dj{k) as the corresponding normalized eigenvector. With this 
notation, e'('^'^~'^j('^)*)(ij(fc) is a standard plane wave solution of p.ip . 

Let us define the matrix D to be the (unitary) matrix with j''th row dj{k). 
This matrix can be used to diagonalize 7i, i.e: 



iuji{k) 







iujj{k) 



D 



(1.2) 



Let e^* denote the propagation operator for (jl.ip . i.e. the operator mapping 
u{x,Q) to u{x,t). In the frequency domain, the propagator can be written as: 



exp 



iuji{k) 







iujj{k) 



iujn{k) 



t D 



(1.3) 



Now consider an initial condition uq{x) localized in frequency about the 
point fco and in position about the point a;o. Consider the solution u{x,t) with 
u{x,t = 0) = ""0(2;). This solution can be written in the frequency domain as: 



ulk,t) — e uo{x)uo{k) 



exp 



iLUi{k) ... 
iujj{k) 

... ILUnik) 




(1.4) 



where u^{k) is the projection of uo{k) onto dj{k) and f{k) denotes the Fourier 
transform of f{x). 

The j-th component corresponds to a superposition of plane waves propa- 
gating with velocity VfcWj(fc). If we consider only regions of fc-space in which 
dkiUjj{k) > 0, then we are considering only waves with a rightward moving 
component. The same can be said about other directions. It is this property 
combined with phase space localization techniques which we will use to filter 
outgoing waves. 
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2 Method 



The methodology of the phase space filter is rather simple. First, given the 
operator 7i, we find the generalized eigenfunctions and the dispersion relation. 
Suppose that the dispersion relations ojj{k) and dj{k) are given. If this is the 
case, then dj{k) is a plane wave propagating with group velocity VkUJ^ik). 



2.1 The Propagation Algorithm 

The propagation algorithm is as follows. We solve (jl.ip on the region [—L — 
w,L + w]^ with any accurate interior solver. At this point, we assume the 
initial condition ^0(2;) is approximately band- limited. We assume that the mass 
of uo{k) for k > fcmax is small. This assumption is common: these are simply 
frequencies too high to resolve with the computational grid. 

We also assume that the frequency content of uq (k) is not localized near the 
regions where the group velocity turns around. I.e., let Z ~ {k : jdk^ujj[k) = 
0} be the set of points where the group velocity vanishes in some direction. Let 
Z be all points a distance fcf, from Z . We assume that for some small J, 

^ ^ 2 ^ 

uo{k) dk<5 (2.1) 

z 

This requirement is necessary due to the Heisenberg uncertainty principle. Very 
roughly, the Heisenberg uncertainty principle says that if the width of the filter 
region is w, then we cannot localize on any region of /c-space smaller than 
(up to constants). To localize only outgoing waves, we would need to isolate 
waves with group velocity pointing outward, while not localizing on waves with 
group velocity pointing inward. Localizing this accurately on waves inside Z 
would be impossible, due to the uncertainty principle. This difficulty can be 
resolved at logarithmic cost, and this is discussed briefly in Section [H 

The propagation algorithm is simple. Fix a time Tstcp < 'fi'/Swinax, with 



sup sup 

3 |fe|<fc„„ 



VkLo,{k) (2.2) 



being the the largest group velocity relevant to the problem. This criterion 
ensures that waves cannot cross the buffer region in a time interval shorter than 
T 

step ■ 

On the time intervals [0,rstcp]7 [Jstop, 2T'stcp]7 ■ • we solve (jl.ip with the 
interior propagator. At times Tstop, 2Tstcp, . . ., we apply outgoing wave filters 
(defined in the next section) to the regions [— L— w, L+w\^\ [— L, L]^ . After 
the application of this filter, all waves which would have reached the boundary 
before a time Tstcp have been removed. Thus, no waves actually reach the 
boundary, and the boundary conditions are irrelevant. 

Algorithm 2.1 TDPSF Propagation Algorithm 

Given an initial condition uo{x), this algorithm calculates Ud{x,t). 
Input: 
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• The dispersion relations and diagonalizing matrices, 0Jj{k) and D. 

• e^*"*, a propagator that accurately solves the interior problem. 

• ^max; the maximal frequency of the problem. 

1. Define T^top = w/ivynax- 

2. Define the approximate solution Ud{x,t) recursively. At times which are 
an integer multiple ofTgtcp, we filter off the outgoing waves: 



Ud{x, (to + l)rstop) 



Ud{x,mTstcp) 



\[{i-o+){i~oz) 

(2.3a) 

The outgoing wave filters are computed using Alaorithm \2.S\. described in 
Section [ 



For other times, we use the given interior propagator: 

Ud{x, TOTstcp + t) = e^''''ud{x, TOTstcp) for r G [0, T^top) (2.3b) 

Ud{x,0) ^ uo{x) (2.3c) 
It now remains to construct the outgoing wave fihers, . 

2.2 Construction of the boundary filter 

Consider a fixed boundary region, say the boundary at xi — L. For a frequency 
fc, if diU}j{k) > 0, then waves with frequency k are outgoing at this boundary. 
The outgoing region of phase space at the right boundary is therefore 

{(f, A?) e X : xi > L and diLUj{k) > 0}. (2.4) 

We now construct a projection onto this region. The Heisenberg uncertainty 
principle makes an exact projection impossible, but we will do the best it allows. 
Extend the box a width w, to be specified shortly. Define the fimction 



N 



xfi^)-[-^) e--^/-^/,(.) (2.5) 

where/,- (a;) = 1 for Xj e [±{L + w/3),±{L + 2w/3)] and Xk G [-L~2w/3,L + 
2w/3] (for k 7^ j), and elsewhere. The parameter a = 0{w/ \n{6~^y^^); a 
precise bound is given in (|2.9p . This ensures that xf{^) < ^ foi' ^ [^L, ±(L+ 
w) or Xk ^ [—L ~ w, L + w]. This function is smooth, and well localized inside 
the buffer region on the j'th sides of the box. 

The set Rj^i — {k G M.^ : dkjOJi{k) > 0} is the set of frequencies with 
the k'th branch of the group velocity pointing right. Due to the Heisenberg 
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uncertainty principle, we cannot project onto this set precisely. However, we 
can approximately project onto most of these wave-vectors. Define 



Rj,LS = {keRk. d{k, Rli) > h} (2.6a) 

The set of vectors within a distance ki, of group velocity dkj oji (x) = is a the set 
of group velocities with motion normal to side j approximately equal to zero. 
The width kh is a buffer to ensure that the frequency spreading caused by our 
spatial localization operators does not cause an error larger than S. Given kh, we 
must also choose a > 0(fcb~^ ln((5~^)^/^) to ensure that the spatial localization 
operators do not spread the frequency content of the solution past the buffer 
region of width kf,. 

In particular, we must choose the buffers widths kb, w and standard deviation 
a to satisfy 



fcb-i (ln(ri)+ln 



^51^5 )) 

w 

<a< , (2.6b) 
" " v/ln((5-i) + iVln(2f77r-i/2) 

which ensures that spreading in frequency does not turn waves around (thus 
minimizing reflection). This is estimated in Section 12.2.11 In particular, note 
that (j2.6bp implies the buffer width w must be at least 0(fcb~^ ln((5~^)) (as 
expected from the Heisenberg uncertainty principle). 
The frequency projection operator is defined as: 



P,-i,^(fc) - ^ e 



2a^^ 



/TT 



-k, a 



... lfl„„,.(fc) 



(2.6c) 



Thus, the operator Pj^^sik) is a smooth projection (in the basis of eigenvectors 
of Ti.) onto wave-vectors propagating rightward. Conjugation of P{k) by D 
moves the projection into the domain of frequency vectors. Finally, we define 
the operator: 

0+ = xti^)D^Pik)Dxt{x) (2.7) 

This operator both localizes in the buffer region at x — L, and projects onto 
waves with group velocity pointing to the right. This operator can be computed 
efficiently and with spectral accuracy as follows: 

Algorithm 2.2 Outgoing Wave Filter Algorithm 

This algorithm applies the outgoing wave filters, i. e. it numerically approx- 
imates (1 — 0^)u{x). 

Input: 

• The dispersion relations and diagonalizing matrices, ujj{k) and D. 



A function il{x). 
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• Indices {j, •) with j G I . . .n and ■ G {+, — }. 

We assume that the operator D^P{k)D is precomputed. 

1. Compute the function xf{x)u{x), and take it's Fast Fourier Transform. 

2. Apply D^P{k)D to the Fast Fourier transform ofxfix)u{x). 

3. Apply the inverse Fast Fourier Transform, and multiply the result by 

4- Subtract the result from u{x), and return the result. 
2.2.1 Choosing the Buffer Widths 

The buffer widths, w and fcf, must be chosen carefully for this algorithm to 
work. There are two competing concerns, namely width of the buffer region 
and frequency spreading, which we must address at this point. 

If w is too small, then the spatial localization functions xfi^) '^iH become 
rougher. But roughness in Xj^i^) "^iH spread the frequencies of the solution 
around. In the frequency domain, multiplication by xfi^) corresponds to con- 
volution, and behaves much like a diffusion operator in the k variable. The 
danger is that the diffusion in k might move mass from the region with positive 
group velocity to the region with negative group velocity. 

However, a larger w will increase the size of the computational box, and 
therefore the computational complexity of the method. Thus, it is desirable to 
take w as small as possible. 

We analyze this as follows. Multiplication by Xji^) corresponds (in the k- 
domain to convolution witlH {wL'^~^ /'i){2^ TT~'^/'^)e-^^'''' . Thus, we can 
approximate the frequency domain operations by: 



I„2_2 , 



The frequency domain operators Pj^i^s{k) are localized on the regions Rj,k,s, 
also with a Gaussian tail of the form (2CT7r~^/^)^e^'"'^°"^ . 

The region Rj^i^s is separated from the region containing incoming waves by 
a buffer of width fc;,. We want to make sure that the spreading in frequency is 
small past this buffer of width k},, i.e.: 



< 5 

\k\=ku 



■^We have suppressed the sinc{...)-factor corresponding to the Fourier transform of Ij{x), 
just for simplicity. This is an overestimate. 
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This can be guaranteed by: 



,,2 tN-1o3N„3N\ \ 1/2 



h-^ [ Inid-^) + ha ( " ^ ^J,, " jj <a (2.8) 

On the other hand, if we take a too large, then the tails of xf{^) '^iU enter 
the computational domain. To ensure that this is minimized, we want to make 
certain that xfi^) ^ ^ foi' ^ € [—L,L]^. By examining the form of (|2.5p . we 
find that this can be accomplished if: 

N 

2 / 2 



or equivalently 



a < , (2.9) 
" v/ln((5-i) +iVln(2a7r-i/2) 



Additionally, to satisfy simultaneously (|2.8p and (|2.9p . we must have that: 

^2_^Ar-l237V^3JVX X 1/2 



-1 (^ln(<5 



fch'^ ln(,5-i)+ln 



^3Af/2 



< , (2.10) 

" v/ln(5-i) + 7Vln(2a7r-i/2) 



This condition demands essentially that w > C'kb~^, where C = 0(ln((5~^)). 
This is precisely what we should expect based on the Heisenberg uncertainty 
principle. 



2.2.2 Powers of 2 

As a practical matter, there is an additional constraint on the buffer width. The 
projections onto outward moving group velocities (the D^P{k)D part of O^) 
are performed using an FFT. The FFT algorithm is fastest when performed 
on a grid of size 2™ in each dimension. For this reason, it is efficient to take 
w = 2™5a; with m [log2 (t«rnm /^a^)! ■ Here, Wmin is minimal w satisfying 
(j2.10p and 5x is the lattice spacing in position. 



2.3 Stability 

The stability of the method is readily proved. The main reason that the algo- 
rithm is stable is simply that it is dissipative: the filtering operator {1 — Of) has 
norm bounded by 1 and can therefore not increase the norm of u{x,t). Addi- 
tionally, all dissipation occurs only at discrete instants of time, which minimizes 
interactions with the propagator. Using these ideas, it is straightforward to 
prove that the algorithm is stable, under the sole assumption that the interior 
propagator is stable. 
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Theorem 1 The Time Dependent Phase Space Filtering propagation algorithm 
is stable if the interior solver is. In particular, we have the estimate: 



\\U4{^,t)\\L^<\\M^)h2 (2.11) 

Proof. The idea of the proof is to show that the numerical solution operator 
at time t can be written as the product of operators of norm 1. Thus energy 
(L^ norm) at time t must be bounded by the initial energy. 
Define the operator U by: 



U 



At time t = rnTgtep +t (with r e [0, Tgtep]), we can write the numerical solution 
Ud{x,t) as: 

Udix,t) = e^'^U'^uoix) 

By the self-adjointness of Hb, =1- If we can show that ||l — 0^\\ < 1, 

then \\U\\ = 1, implying: 

\\ut,ix,t)\\L^ < lle^^^ll \\U\r\\uo{x)\\^^ <l-l^-\\Mx)h^ = \\Mx)\\l^ 

and stability is proved. 

We first show that (T{Of) C [0,1]. Recall that Of is defined as Of = 

xf(x)D'<P{k)Dxf{x). Note that D'<P{k)D is diagonalizablc (to P{k)) and 
each diagonal entry is contained in [0, 1] for each k. Therefore, D'^P{k)D 



< 1 



and a{D''P{k)D) = a{P{k)) G [0, 1]. This implies that < 1. Now write: 



{f\Off) = {xt{:r,)f{x)\D'^P{k)Dxf{x)f{x 

'Dxf{x)f{x)\P{k)Dxfix)f{x)} = (g\P{k)g 



where g = Dxj{x)f{x). Since P{k) is a positive matrix for each fc, positivity 
follows. Since is a positive operator with norm bounded by 1, has 
spectrum in [0, 1]. The spectral mapping theorem implies that a{l — O^) C 
1 — [0, 1] = [0, 1], implying that ||l — O^H < 1, and stability is proved. □ 



2.4 Tangential Waves 

The filter described in this section is not optimal. In particular, some outgoing 
waves that should be filtered are not, primarily waves which are outgoing, but 
nearly tangentially to the boundary. 

While it is impossible to completely resolve this issue (due to the uncer- 
tainty principle), one can improve the situation by using a better phase space 
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localization scheme. One can build phase space projections using framelets (lo- 
calized wavepackets) fine tuned to the problem of interest. For the Schrodinger 
equation, canonical coherent states (the Gaussian Windowed Fourier transform) 
frame [53] is a natural choice. Other equations require different frames (e.g. 
curvelets or wave atoms for hyperbolic waves 6J). 

The cost of this approach is extra programmer time and (usually) an increase 
in computational complexity by a constant factor. 



3 Examples 

3.1 A warm-up: 1-dimensional Schrodinger equation 

Phase space filters were originally developed for the Schrodinger equation {u{x, t) 
is a scalar field, and Ti, = lA). For the Schrodinger equation, the phase space 
filter takes a particularly simple form. There is no diagonalizing operator, and 
the group velocity is simply k. Thus, outgoing waves at the right boundary are 
simply waves of positive frequency (negative frequency on the left boundary). 
In this case, P{k) takes the simple form: 

Thus, the filtering operator is simply 

(and similarly on the left side). For simplicity, we take fc;, = 0. 

We solved this example on a lattice of 1024 pts, with 5x = 0.1 and initial 
condition 

u{x,t = 0) = ^e-" /^-^ 

We measured the errors for various values of 6, taking a = 1. The results 
are plotted in Figure [T] Note further reducing 5 does not reduce the error 
significantly beyond 10~*; we believe this is due to floating point errors. 



3.2 Euler equations, linearized around a jet 

We now move on to a more interesting example. 

We consider now the two-dimensional Euler equations linearized around a 
homogeneous jet flow in the xi direction. The vector u{x,t) is 3 dimensional, 
with ui representing the pressure change, W2,3 representing the velocity field in 
the xi^2 directions (respectively). In this case, Ti. takes the form: 



n = 



Md^^ 
Md,, 



(3.1) 
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Figure 1: The relative error for the 1-dimensional test of the Schrodinger equa- 
tion, as a function of k and the parameter 6. 



Here, < M < 1 is the mach number. Ti, has eigenvalues uji{k) — Mki 
uJ2{k) = Mki — \k\ and L03{k) — Mki. The diagonalizing matrix is 



D 



V2\k\ 



~\k\ ki k2 
\k\ ki k2 
-%/2fci V2k2 



(3.2) 



When M 0, this is the classical example of an equation for which the PML 
becomes unstable 

We solved the Euler equations on the region [—32,32]^ with lattice spacing 
Sx = 0.125 (for a total of 512^ lattice points). Such a setup is valid for spatial 
frequencies up to fcmax = tt/Sx — 25.1. We solved the system using the Fast 
Fourier transform to compute (jl.3p . which is accurate to machine precision 
provided no waves reach the boundary ( |22l Theorem 4.1] or [3T1 Theorem 
A.l]). This accuracy is independent of time-step, which we took to be 6t = 0.25 
in order to get a watchable frame-rate in the plots. The phase space filter region 
was taken to have width 16 (128 lattice points), with Tgtop = 1-5, and a = 1.0. 
The initial condition was taken to be 



with 



ui(a;,t = 0) = r^e^''^/^ cos{Kr) 
r = V(a;-8)2-hy2. 



(3.3a) 
(3.3b) 
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Figure 2: The pressure field at various times. The arrows indicate the di- 
rection of the velocity field. Note that the phase space filter was applied be- 
tween t = 9.75 and t = 10.0. This simulation is available in movie form from 
http://cims.nyu.edu/ ^stucchio/software/kitty/ demos/euler_tdpsf.mpg[ 
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Figure 3: The relative error (measured in various norms) as a function of the 
frequency of the initial condition. 



with the K varying from 1 to 20 (the frequency range of the problem). This 
yields an initial condition with frequency localized near |fc| = K. The results 
were compared to another simulation on a large box. The result is that for 
K > A, error of lO^'^ (relative to the initial condition) is achieved up to 
t = 50 (see Figure [3]). The error at low frequencies can be dealt with by various 
means, see Section |4] and [21] . 



3.3 Maxwell's Equations in an Orthotropic Medium 

Let E be the electric field and H = (xB with B the magnetic field. Let e and ^ 
be the electrical and magnetic permeability's of a medium. 

In order to bring Maxwell's equations to symmetric form, we introduce the 
auxiliary variable u = {^JJiH , y/eEy . With this formulation, the Hamiltonian 
for Maxwell's equations can be written in block form: 



n 







-1/2V 



X ^ 



-1/2 



-Ai-i/2V X e 




-1/2 



where Vx is interpreted as a matrix. 



(3.4) 



Vx = 





9. 



-d, d, 





-dy dx 



y 




The symmetry of e, ^ and Vx implies that Ti is skew-adjoint. 
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To further simplify the system, we make the following additional assump- 
tions. We assume fi = 1. We assume the system is z-independent, i.e. dz = 0, 
and we can restrict ourselves to two-dimensional simulations. Lastly, we simplify 
the electrical permittivity: 



1 b 
b 1 
0c 



This is the simplest possible birefringent system. With the variables 

/ = (l/2)(VTT6+x/T~5) 
9 



(3.5) 



= (i/2)(-VTT6 + VT^), 



we can write the dispersion relation as: 



0. 



(3.6a) 
(3.6b) 
(3.6c) 



with 



D = 



V2|fe| 
k2 



\/2|fe| 






fcl 








fc2 









-2-1/2 

2-1/2 






fk2-gki 

V2E{k) 
fk2-gki 






/fcl 



-gk2 



V2E(k) 





/fcl 



-gfc2 



E{k) 



V2E(k) 
/fci-gfc2 
V2E(k) 



/fc2-gfci 
E{k) 



2-1/2 
2-1/2 









(3.7) 



with E{k) = ^{p+g^){kl + kl)-Afgk^k2. 

This can be further simplified by noting that ui,U2 and uq (corresponding 
to Bx,By,Ez are uncoupled to U3,U4 and M5. The ui,U2 and uq modes (cor- 
responding to Bx,By and Ez) are actually isotropic waves, so we will ignore 
them. We therefore restrict consideration to U3,U4 and U5 (1*3 corresponds to 
Bz, and M4_5 correspond to to linear combinations of E^ and Ey). 

We solved the reduced system of Maxwell's equations with the same param- 
eters as in Section 13.21 with the same initial condition (replacing ui by in 
(|3.3ap ). The anisotropy parameter b was chosen to be 0.25. The results are 
comparable to those for the Euler equation, see Figures |4] and [31 



3.4 Long Time Stability 

In Theorem [T] of Section 12.31 it was proved that the phase space filtering algo- 
rithm is stable. To demonstrate the validity of the theorem, we ran a simulation 
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Figure 4: The magnetic component Bz of the electromagnetic pulse at var- 
ious times. Note that the phase space filter was applied between t = 
9.75 and t — 10.0. The non-radial shape of the wave is due to the 
anisotropy of the medium. This simulation is available in movie form from 
|http://cims. nyu.edu/~stucchio/software/kitty/demos/maxwell-tdpsf.mpg; 
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Errors for Maxwell's Equations 




Frequency 



Figure 5: The relative error (measured in various norms) as a function of the 
frequency of the initial condition. 



of the Euler equations and Maxwell's equations up to time t — 2000. While we 
can not determine the accuracy over such long time intervals (the reference sim- 
ulation would require an extremely large box) , we can study the growth of the 
L'^-norm. 

The results of such a simulation are plotted in Figure [Sj They indicate that 
Theorem[T]is correct, and that the mass of the solution decreases monotonically 
with time. 



4 The Low Frequency Problem 

As is apparent from Figures [3] and O the phase space filtering approach does not 
work well for waves with low frequency. The reason for this is that to localize in 
frequency, the region in which one works must be 0(1) wavelengths long. The 
simplest remedy is to increase the width of the filter. If the smallest frequency 
relevant to the problem is kb, then the width of the buffer is 0{kb~^), which 
means that the computational cost is of order 0{kb~^). 

This problem can be remedied by a somewhat more involved method, which 
has been implemented for the Schrodinger equation [21j . The essential idea is 
to increase the width of the box, but reduce the sampling rate on the extended 
region. Then high frequency waves are filtered at the edge of the highly sampled 
region, and low frequency waves are filtered on the edge of the coarsely sampled 
region, but using a wider filter capable of resolving low frequency waves. 
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Long Time Stability 




2000 



Time 



Figure 6: The norm of solutions of the Euler and Maxwell equations as a 
function of time. In both cases, we took K — IQ (recall (|3.3aP ). The Mach 
number was M = 0.5 for the Euler equations and the anisotropy was h = 0.25 
for Maxwell's equations. 
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With this method, even though the box has width 0(fcf,~^), the number of 
samples required is only 0(log(fcmax/fcb)) (computation time scales similarly, 
up to logarithmic prefactors). This aUows resolution of outgoing waves at low 
frequencies in logarithmic rather than linear cost. It is also argued heuristically 
in [5T] that this computational complexity is close to the best possible. 

Acknov^rledgements: We thank P. Petropoulos for pointing out to us the 
problem of PML instability, and R. Goodman for a useful discussion on lattice 
waves. 
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